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Abstract: In ghost imaging schemes information about an object is 
extracted by measuring the correlation between a beam that passed the 
object and a reference beam. We present a spatial averaging technique that 
substantially improves the imaging bandwidth of such schemes, which 
implies that information about high-frequency Fourier components can 
be observed in the reconstructed diffraction pattern. In the many-photon 
regime the averaging can be done in parallel and we show that this leads 
to a much faster convergence of the correlations. We also consider the 
reconstruction of the object image, and discuss the differences between a 
pixel-like detector and a bucket detector in the object arm. Finally, it is 
shown how to non-locally make spatial filtering of a reconstructed image. 
The results are presented using entangled beams created by parametric 
down-conversion, but they are general and can be extended also to the 
important case of using classically correlated thermal-like beams. 
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OCIS codes: (270.0270) Quantum optics, (110.0110) Imaging systems, (190.0190) Nonlinear 
optics, (100.0100) Image processing. 
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1. Introduction 

Lately a substantial effort has been put into understanding the physics behind "ghost" imaging 
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], also called entangled imaging, two-photon imaging, 
coincidence imaging, or correlated imaging. The ghost imaging schemes are based on two 
correlated beams typically originating from parametric down conversion (PDC). One beam 
travels a path (the test arm) containing an unknown object, while the other beam is sent through 
a reference optical system (the reference arm). Information about the object is then obtained 
by measuring the spatial correlation between the beams. This two-arm configuration allows for 
obtaining different kinds of information by solely operating on the reference arm optical setup 
while keeping the test arm fixed. In particular, both the image (near-field distribution) and the 
diffraction pattern (far-field distribution) of the object can be measured [9]. 

The early studies concentrated on the working principles of the ghost imaging scheme, both 
in terms of basic formalism [1, 2] as well as the experimental implementation [3, 4, 5]. A 
general theory has been developed that discusses how to extract the information, not only in 
the coincidence counting regime [6, 7, 8], but also in the high gain regime [9, 10, 11, 12] 
where recent experimental results showed nonclassical spatial correlations for for high-gain 
PDC [14], which promises well for an efficient quantum ghost imaging protocol. Most recent 
discussions have instead focused on whether entanglement is necessary for extracting the infor- 
mation [6, 7, 8, 9, 10, 11, 13] and references therein. The present paper focuses on how to op- 
timize any source, entangled or not, to give as much image information as possible. The issues 
of imaging bandwidth and image resolution have been taken up previously [9, 10, 11, 12, 13]. 
Even in an ideal ghost imaging system the imaging bandwidth is limited by the source generat- 
ing the correlated beams, characterized by a finite extension of the spatial gain in Fourier space 
^source: beyond this value the Fourier frequencies are cut off and this information is lost in the 
reconstruction of the object diffraction pattern. The image resolution is instead limited by the 
coherence length jc co h, which is given by the characteristic width of the near-field correlation 
function. 

The question is: can we increase the source cutoff value? In the case of PDC the cutoff 
is determined by the phase-matching conditions [15, 10, 12] while in the case of a classically 
correlated field generated from chaotic thermal radiation, the bandwidth is roughly given by the 
inverse of speckle size found from the self-interference of the near field [16] (see also [10, 13]). 
Thus, the cutoff is an inherent property of the source. However, in this paper we discuss a spatial 
averaging technique which circumvents this cutoff and improves the imaging bandwidth of the 
system, regardless of how the correlated beams are created. The spatial average is performed 



over the position of a pixel-like detector located in the test arm after the object, exploiting that 
each position of this detector gives access to a different part of the diffraction pattern through 
the correlations. Thus, making an average over all possible positions of the test detector it is 
possible to substantially extend the imaging bandwidth of the scheme. 

The spatial averaging technique works particularly well in the high-gain regime where many 
photon pairs are generated in each pump pulse. Thus the test arm has many photons per pulse 
transmitted by the object, (in contrast in the low-gain coincidence counting regime only one 
photon at a time is impinging on the object, and either it is transmitted or it is not), and therefore 
at the measurement plane they are scattered over the entire transverse plane. The information 
about the object is then extracted by spatially correlating the intensities recorded in the test 
and reference arm. In the low-gain regime this implies registering coincidences while in the 
high gain regime successive sampling over repeated shots of the pump pulse is used. Since the 
spatial averaging technique employs an average over the position of the test detector, within a 
single shot we can get information about the image from all the test detector positions in the 
transverse plane containing photons. This implies that in the high gain regime a much faster 
convergence rate is obtained using the spatial averaging technique. 

The spatial averaging technique was already introduced in Ref. [12] in the case where homo- 
dyne detection was used. Since the homodyne measurements give access to orthogonal quadra- 
tures, an increased image resolution can be obtained by using the spatial averaging technique to 
measure the diffraction pattern and then use an inverse Fourier transform to obtain the image. 
In the present work we only have access to the field intensity, so we cannot use this method 
to get an increased image resolution. However, we may use a bucket detector in the test arm 
when we want to observe the image, which was not possible in the homodyne scheme due to 
phase-control problems, and it turns out to make the imaging system incoherent. We discuss 
the benefits and disadvantages in this respect. We will also extend the discussion of [12] and 
give a more detailed and general picture as to why the spatial averaging technique works, and 
as to how much the convergence rate can be improved. 

The majority of the results of this paper hold for ghost imaging schemes in general, however 
we use in the following a model for PDC as the source for the correlated beams. 

The paper starts by presenting the analytical results in Sec. 2, and the spatial averaging 
scheme is introduced and discussed. The numerical results are contained in Sec. 3, going be- 
yond the approximations made in the analytical section and validating the results by using very 
realistic parameters of current experiments. The paper is summarized in Sec. 4. 

2. The model and analytical results 

We consider the PDC model for the three-wave quantum interaction inside a %^ nonlin- 
ear crystal previously discussed in detail in Refs. [12, 15]. The crystal of length / c is cut 
for type II phase matching, and the model consists of a set of operator equations describ- 
ing the evolution inside the crystal of the quantum mechanical boson operators aj(x,z,t) for 
the signal (J = 1) and idler (J = 2) fields, obeying at a given z the commutator relations 
[ai(z,x,f ),aj(z,x',?')] = 5,-/5(x — x')5(f — f'), i,j = 1,2. In the stationary and plane-wave pump 
approximation (SPWPA) unitary input-output transformations can be written relating the field 
operators in q and Q. space at the output face of the crystal bj(q,£l) = aj(z = / c ,q,£2) to those 
at the input face a™(q,i2) = aj(z = 0,q,£2) as follows 

b j (q,Q.) = U j (q,Q.)a i ] 1 (q,ty+Vj(q,Cl)af(-q,-Cl), j^k=l,2- (D 

The gain functions Uj and Vj can for example be found in [12, 15]. 

The system setup we consider is the same as in Refs. [9, 10, 1 1, 12] and has the characteris- 
tic two-arm configuration of the ghost imaging scheme, see Fig. 1. The signal and idler exiting 
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Fig. 1. System setup. In (a) the setup for reconstructing the diffraction pattern of the object 
is shown, where the reference arm is in the f-f setup. In (b) the reference arm is changed to 
a telescope setup, which is used to reconstruct the object near field. PBS: polarizing beam 
splitter. L: lens of focal length /. T^: object. D\, Di'- arrays of detectors. 



from the crystal are separated with a polarizing beam splitter (PBS). The object is located in the 
test (signal) arm immediately after the PBS, and the transmitted field then propagates through 
an f-f lens system: a lens placed distance / from the object as well as to the measurement plane, 
where D\ (an array of pixels or a CCD camera) is recording the intensity of the field. Thus D\ 
is in the focal plane of the lens, but it is important to stress that it cannot on its own measure 
the interference pattern of the object because the signal beam itself is incoherent. This test arm 
setup is kept fixed, while in the reference (idler) arm two different configurations are used. In 
Fig. 1(a) the reference arm contains also an f-f lens system, and by recording the idler field 
intensity with (another CCD camera) and correlating it with the signal intensity, informa- 
tion about the object diffraction pattern can be extracted from the correlations. Conversely, in 
Fig. 1(b) a so-called telescope setup is used consisting of two lenses with focal length / placed 
4/ from each other and 2/ from the PBS and D2, respectively. This setup allows for retrieval 
of the object image distribution from the correlations. 

The analytical results for the correlations in the two cases of Fig. 1 were derived in 
Ref. [9, 10], and we review them in the following. For simplicity we neglect for the moment the 
temporal argument, which corresponds to using a narrow frequency filter after the crystal. The 
propagation of the fields from the PBS to the measurement planes are described by the kernels 
hj(xj,x'j). The fields at the measurement plane Cj(xj) are then found by the Fresnel transfor- 
mation cj(xj) = f dXjhj(xj ,x'j)bj(x'j) +Lj(xj). Ljixj) represent the losses in the propagation 
which are proportional to the vacuum field operators and are uncorrected from bi{Xj). The field 
intensities at the measurement planes Ij(xj) = Cj(xj)cj(xj) are detected with Dj, and result in 

the correlation (/] (xr )/2(x2)) = (c}(xr )ci(xi)c2(x2)c2(x2)). The information about the object 
is obtained by subtracting the background term to obtain the intensity fluctuation correlations 

G(x 1 ,x 2 ) = (/ 1 (x 1 )/ 2 ( X2 )>-(/ 1 (x 1 ))(/ 2 (x 2 )). (2) 

It is then straightforward to show the following essential result [9, 10] 



G(xi,x 2 ) 



dx[ / dx^ 1 (x 1 ,x 1 ')/ l2 (x2,x^)r(x;-x^n = 0) 



(3) 



Here r(xi,X2,n) is the near-field correlation at the crystal exit, and in the SPWPA and with a" 



in the vacuum state it is found from Eq. (1) as 



r(x 1; x 2 ,£2) ee J dze-^(b l (x,,t,)b 2 (x 2 ,t 1 +T)) = I t^etoto-^faG), (4) 

where we have introduced the gain function y(q,£2) = U\ (q,£2)V 2 (— q, — £2). Since Eq. (4) 
is a function of xi — x 2 (because of the SPWPA) we will use the notation r(xj — x 2 ,£2) = 
r(xj , X2 , £2) in the following. 

We should mention that when temporal argument is taken into account, Eq. (3) is no 
longer valid. We consider the intensities averaged over the detection time Td as = 
Tp l f T 6tCj(Xj,t)cj(Xj,t). When T^, is much larger than the coherence time of the source T co h 
(which for PDC is typically the case) the following expression holds instead [11] 



G(Xl ' X2)= 7bi2^ 



d X ; / dx^ l (x l ,x;)/ 12 (x 2 ,x 2 )r(x;-x^ ) £2) 



(5) 



This form will be used later for comparing with the numerics. 

We fix the test arm once and for all as shown in Fig. 1, so hi(xi,x[) <x e _,Xl ' x i A/ '^7' bj(xJ), 
where k is the free space wave number of the degenerate down-converted fields. To extract 
information about the diffraction pattern the reference arm is set in the f-f setup of Fig. 1(a), 
h2(x2ix' 2 ) <x e~ , * 2 '*2 lc /f, and Eq. (3) then straightforwardly implies the correlation [9, 10, 11] 

G f (x l ,x 2 )°c\ r (-x 2 k/f,£l = 0)f obj [(x 1 +x 2 )k/f}\ 2 , (6) 

where 7obj(q) = / ||e _ ' q ' x Tobj(x) is the amplitude of the object diffraction pattern. The sub- 
script "f" denotes that the reference arm is in the f-f configuration. The correlation provides 
information about the diffraction pattern of the object if we fix X] and scan x 2 , but since the 
gain also depends on x 2 there is a limit to the information we can extract. Precisely, the gain 
y(— x 2 £//,£2 = 0) introduces a cutoff of the reproduced spatial Fourier frequencies at a certain 
characteristic value which we denote dqpoc', the imaging bandwidth of the PDC source. 

We will now show how to circumvent this source-related limitation on the imaging band- 
width. As previously shown [12] we may in a suitable way average over the position of the test 
arm detector Xi: First, a change of variables is introduced as x = xi +x 2 , and then an average 
over xi is performed. We then obtain 

G fiSA (x)EE /dxiG f (xi,x-Xi)« /d Xl | 7 [(xi-x)A://,£2 = 0]f ob j(xV/)| 2 



= \T oh] (xk/f)\ z J dx, |y[( Xl -x)k/f,Cl = 0]\ £ ~ const x |r obj (xA://)| 2 (7) 

where the subscript "SA" indicates that a spatial average has been carried out. The final approx- 
imation in (7) is that | y(q, £2 = 0) | 2 is a bound function inside the detection range of xi implying 
that the integral evaluates to a constant. Thus, there is now no gain cutoff of the diffraction pat- 
tern, so the imaging bandwidth becomes practically infinite (apart from limitations arising from 
the finite size of the optical elements). 1 Note that this average over xi does not correspond to 
D{ being a bucket detector. Instead, the change in variables x = xi +x 2 implies that the spa- 
tial averaging technique is a convolution of the signal and idler intensity fluctuations, which in 
practice is easily calculated using fast Fourier transform techniques. 

1 We should mention that when x is taken to the boundary of the detection range, the integral is no longer a constant. 
Thus at the boundaries the bandwidth slowly dies out, but it is a quite small effect only affecting a range of Sqmc 
there. In our numerical simulations we use periodic boundary conditions, so this limit does not come into play. 
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Fig. 2. The spatial averaging technique explained by the Klyshko picture. The lower sketch 
shows the unfolded version of Fig. 1(a). The red line represents the ray from the detector at 
position Xj as it travels through the system and reaches its conjugate point in the reference 
plane on the right side. The upper plot shows the source gain curve in Fourier space and 
the position of a typical diffraction pattern for the current value of jc9 ; only the black part 
is amplified while the gray part is not amplified. As the test detector position is changed 
to x'j + Ax the dashed red ray comes into play, and the diffraction pattern in the upper 
plot moves a corresponding amount (also dashed). The gain instead remains the same, and 
hence a different part of the diffraction pattern is amplified. (Movie size 449 KB). 



In Fig. 2 we give an intuitive explanation of why the bandwidth is extended by averaging over 

Xi, and it is based on the Klyshko picture [1] (see also [17]). This states that the propagation 

in the two distinct arms can be viewed in an "unfolded" scheme having the nonlinear crystal 

as the hinge point of the signal and idler rays, i.e. the crystal can be treated as a geometrical 

reflection mirror. This is a consequence of momentum conservation in the PDC process, so a 

signal ray with momentum q' implies an idler ray with momentum — q'. Thus, we may treat 

the scheme starting from the signal detector with the signal ray passing backwards in time 

through the test arm, until it encounters the crystal and becomes the idler ray traveling forward 

in time. The point is that the action of the crystal is to fix the idler ray propagation direction 

through the momentum conservation, and in the unfolded version of Fig. 2 the idler ray is 

a straight-line continuation of the signal ray. According to the Klyshko picture the detector 

in the test arm acts as a point-like source emitting a spherical wave from Xi. This wave is 

transformed by the f-f lens system into a tilted plane-wave with transverse wave vector q = 

—X\k/f which then after impinging on the object is diffracted. Inside the crystal the diffracted 

signal ray is "converted" into the idler ray with "transmission" function y(q, D. = 0) determined 

by the phase-matching conditions. This transmission function y(q,£2 = 0) is more precisely 

amplifying certain components of the signal ray and damping others (gray area in Fig. 2), and 

this is the origin of the infamous cutoff. Therefore, as the idler ray now propagates forward 

through the f-f lens system in the reference arm we can with D2 only detect a bandwidth- 

1 ~ 1 2 

limited version of the diffraction pattern y(— X2k/f,Q. = 0)r o bj [(X2 +xi)k/f] \ . The spatial 

averaging technique exploits that as we move to a new detector position x j + Ax the center of 

the object changes correspondingly to X2 = — x^ — Ax, while the gain remains fixed [due to the 

structure of Eq. (6)]. Thus, by changing the test detector position, we measure a different part 



of the spatial Fourier spectrum, as shown in Fig. 2. Therefore, by performing a suitable average 
as described in Eq. (7) we cover the entire Fourier plane, and effectively the bandwidth has 
become unlimited. 

Additionally, the spatial averaging technique gives an increased convergence rate of the cor- 
relation. This is again related to the fact that when the test detector position is changed from 
Xj to Xj + Ax the gain does not change position but the diffraction pattern does. Assuming that 
the PDC gain extension is much larger than the extension of a pixel, the shifted diffraction 
pattern at Xj + Ax overlaps quite substantially with the previous one. Thus, as xi is scanned 
a given position of the diffraction pattern has as many independent contributions as there are 
independent modes inside the gain bandwidth, which as a good estimate is given by the ratio of 
the PDC bandwidth dqpoc with the pump bandwidth 8q p = 2/wq [15], where wq is the pump 
waist. A measure of the speedup in the correlation convergence in each transverse dimension is 
therefore given by 

Psa = Sq PDC /8q p . (8) 

In the simulations shown later concerning the convergence rate, we used a temporal grid that 
corresponds to a 22 nm interference filter. In that case Sqpoc — 6go> where qo = \fk/T c is tne 
natural bandwidth of PDC at a given temporal frequency [15]. The pump size in the numerics 
was chosen so 8q p — go/ 18 (corresponding to a pump size of 600 |im), implying we should 
expect a convergence rate speedup of Psa — 10 2 . 

Note that instead of fixing xi and scanning X2, we may scan Xi and fix X2. In this case Eq. (6) 
reveals that the gain no longer limits the imaging bandwidth, and therefore it is in principle 
not necessary to do a spatial average to overcome the gain limitation. The physical explanation 
behind this is again found from the Klyshko picture: fixing X2 at a suitable position (i.e. at the 
position of maximum gain, cf. Fig. 2), scanning xi will move the diffraction pattern seen at this 
position over the whole range giving an unlimited imaging bandwidth. However, scanning Xi 
and fixing x? does not benefit of the improved convergence. To achieve this a spatial average 
should be done, whereby the method amounts to the same operation as described previously. 

Let us now turn to the problem of reconstructing the object image. Keeping the test arm fixed 
and changing the reference arm to the telescope setup [see Fig. 1(b)], /i2(x2,X2) = 5(xi — Xj) 
and Eq. (3) becomes [9, 10, 1 1] 



Gt(xi,X 2 ) oc 



dxjr( x ; - X2 ,n = o)r obj (x() 



2 



(9) 



|r obj (x 2 )| 2 | 7 ( Xl /://,£2 = 0)| 2 . (10) 



Thus, fixing Xi the object can be observed by scanning X2. The approximation that leads to 
Eq. (10) holds as long as the smallest length scale of the object is larger than the coherence 
length x co h of T(x[ — X2,fi = 0). This is because r(x( — X2,i2 = 0) is nonzero in a region 
of the size x co h around xj = X2 and thus r o bj changes slowly with respect to this function 
and may consequently be pulled out of the integration. In general, Eq. (9) is a convolution 
between the correlation function and the object, and hence the width x co h of the near-field 
correlation function r(x[ — X2,£2 = 0) determines the resolution of the reconstructed image, 
and has typically a value of x co h = 1/qo- 

In reconstructing the image a technique corresponding to the spatial average done for recon- 
structing the diffraction pattern would result in a largely increased image resolution. Unfortu- 
nately, it is not possible to carry out such a spatial average to achieve this. However, if we make 
a simple sum over all the positions of D\ (corresponding to D\ being a bucket detector), instead 



of Eq. (9) we have the following exact expression 2 



G T (x 2 ) 



/ 



dxiG T (xi,x 2 ) 



/ 



dx 1 |r(x 1 -x 2j n = o)| 2 |r obj (x 1 ) 



2 



(ID 



where the bar denotes that D\ is a bucket detector and therefore that xi has been integrated 
out. Eq. (11) has the form of an incoherent imaging scheme, while Eq. (9) has the form of a 
coherent imaging scheme (the same conclusion - that using a bucket detector can make the 
imaging incoherent - was reached in Ref. [8] in the low-gain regime). As we shall see later, the 
incoherent form has both advantages and drawbacks compared to the coherent case. 

3. Numerical examples 

The imaging performance of the system was checked through numerical simulations of the 
model described in detail in Ref. [15]. It includes spatial and temporal dispersion up to second 
order, as well as a Gaussian shape in space and time of the pump beam. The simulations with 
1 transverse dimension (ID) were calculated including the temporal argument (using a grid of 
N x = 512 and N t = 64), while the ones with 2 transverse dimensions (2D) had a more qualitative 
nature since they neglected the time argument (a grid of N x = iVy = 256 was used). 3 Each 
pump shot was propagated along the crystal in N z = 200 steps. The Gaussian pump profile 
had a waist wq = 600 jxm and a duration of To = 1.5 ps, and the other parameters were as 
in [12, 15] chosen to closely mimic that of an l c = 4 mm long BBO crystal used in a current 
experiment in Como [ 14] . The characteristic space and time units of the PDC source are x con = 
1/qo — 16 jUm and T co h = 1/£Iq — 0.96 ps. The pulsed character of the pump is important 
because typically the CCD cameras used in experiments have very slow response times. If 
the measurement time becomes too long with respect to T co h, the visibility of the correlation 
becomes very poor (see also Refs. [6, 10, 1 1]). This only affects the case when the intensity is 
detected in the high gain regime, i.e. when the background term in Eq. (2) is substantial, while it 
does not affect the coincidence counting regime or the case when homodyne measurements are 
performed as in Ref. [12]. Note that in the telescope setup the performance was optimized by 
taking the imaging plane of the telescope setup inside the crystal by the amount Az = — (l/«i + 
1 /«2)tanh((7pZ c )/c7p [12, 15], where tij are the refractive indices, and <7 p is a gain parameter. In 
this way, the quadratic phase dependence of the gain is cancelled. For more technical details on 
this and the numerics we refer to Refs. [12, 15]. 

Let us show quantitatively how the object diffraction pattern is reconstructed by using a 
double slit as an object in the setup of Fig. 1(a). The correlations as calculated in ID both 
with and without the spatial average are presented in Fig. 3(a)-(b). The correlation without 
spatial average in (a) clearly suffers from a limited bandwidth, and from the animation a slow 
convergence rate is evident. In contrast, the correlation with spatial average in (b) is able to 
reproduce the entire spectrum of the diffraction pattern and after much fewer repeated pump 
shots. To check the convergence rate, we have calculated the root-mean-squared error of the 
correlation G" (after averaging over n number of shots) relative to the analytically calculated 
correlation Gf. Gf was calculated by using semi-analytical methods including the analytical 
SPWPA gain as well as integrating over time [see Eq. (5)]. The error is then evaluated as 



2 The expression is exact because no approximations have been made about object length scales vs. jr co h- 
3 This was done to save CPU time since many thousands of repeated pump shots were needed for the correlations 
to converge, and corresponds to using a narrow bandwidth interference filter after the crystal. In fact, the temporal grid 
acts just as an interference filter by providing a cutoff in frequency space. In the simulations with time the cutoff was 
chosen so it corresponds to having a 22 run interference filter after the BBO crystal, while the simulations neglecting 
time are equivalent to placing a < 1 nm interference filter after the crystal. 
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Fig. 3. Reconstructing the diffraction pattern using the f-f setup in the reference arm. (a) 
shows the correlation from a fixed x\ after averaging over 10000 shots while (b) shows the 
correlation using the spatial averaging technique and averaging over 1000 shots. In (a) and 
(b) black curves are numerics, while dashed red curves are . In (a) the blue curve is 
the analytically calculated correlation in the SPWPA in arbitrary units [Eq. (6) including 
integration over time, as shown in Eq. (5)]. (c) shows e(n) from Eq. (12) in a double- 
logarithmic plot using the spatial averaging technique (full) and a fixed detector (dashed), 
while the red curves are fits to Eq. (13). Double slit parameters: aperture 14 /im (5 pixels) 
and aperture center distance 87 flm (30 pixels). [Movie sizes: (a) 1 MB and (b) 933 KB]. 



where at each shot the maximum of G" is rescaled to the maximum of Gf. Figure 3(c) shows 
e (n) and clearly the fixed detector case (dashed) converges much slower than the spatial average 
(full). The red curves are fits to the function 

e fit («) = 1/2 + (13) 

with n being the number of shots. The factor do has the dimension shot - 1 and the ratio between 
afo's obtained using the spatial average and using fixed detector should then give an indication 
of the increased convergence rate. From the fits of Fig. 3(c), cIq AAD / dQ Ked ' lD ~ 125, which cor- 
responds well to the value Psa — 10 2 estimated in the previous section. The factor d\ indicates 
the error remaining after convergence, but in this case the correlations are not converged suf- 
ficiently within the number of shots shown. Now, according to what was predicted in Eq. (8), 
using the spatial average in 2D should benefit with a factor of Psa in convergence speedup com- 
pared to using the spatial average in the ID case (due to the extra dimension to average over). In 
Fig. 3(c) the result of such a 2D simulation is also shown, 4 and the improved convergence rate 

4 This was a full 3+1D simulation with N x = 256, N v = 128, N. = 200 and N, = 32. The slits had the same extension 
in the x-direction as in the ID case and were 580 /im (90 pixels) long in the y direction. 



of the 2D simulation is evident: from the fits d Q ' /d ' ~ 73, again corresponding well to 
the predicted estimate of Psa — 10 2 . A minor problem in this comparison is that the 2D results 
saturate very quickly to a rather high error level of 1%, while the ID results go as low as 0. 1%. 
This difference turned out to be caused by the Gaussian shape of the pump field and the object 
extension; the object is very localized in the x direction and therefore we get a low error rate 
in the ID case. However, in the 2D case the object is quite extended in the y direction and the 
average error reported in the £ value is therefore higher. This quick saturation in 2D gave some 
numerical problems in fitting well the data to Eq. (13), so the do value obtained should be taken 
cautiously. However, there is no doubt about the main point: going from ID to 2D the spatial 
averaging technique speeds further up. 
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Fig. 4. Reconstructing the diffraction pattern of the object shown in (a) which has spatial 
Fourier components both inside and outside the PDC bandwidth. The real part of the inverse 
Fourier transform of the reconstructed diffraction pattern using the f-f setup is shown with 
(b) fixed xi and (c) using the spatial averaging technique. 35000 shots were used. 

Let us show an example where the reconstructed diffraction pattern using the f-f setup leads 
to completely different results depending on whether the spatial average is applied or not. The 
object is shown in Fig. 4(a) and is a 2D mask r o bj(x) = [1 +cos(xgo)][l + cos(3y^o)]/4 (con- 
sisting of 4 main sidebands, 2 located at q/qo = ±e v and 2 located at q/qo = ±3e v ). We have 
chosen the peaks in the y-direction so they lie outside the imaging bandwidth of the source, 
while the peaks in the x-direction lie inside the bandwidth. In the numerics the diffraction 
pattern was reconstructed from the correlation with xj fixed as well as employing the spatial 
averaging technique. Instead of showing these data, we have in Fig. 4(b) and (c) shown the 
results after taking the inverse Fourier transform (IFT) of the reconstructed diffraction patterns. 
In this way we can reconstruct to some extend the nature of the near field mask (some phase 
information is lost, obviously, but in this simple case the lost phase information is not so im- 
portant). Figure 4(b) shows Re(IFT[Gf(xi,X2)]), i.e. without spatial average, and only the roll 
pattern in the .^-direction is reproduced. Instead, the correct object is reproduced by the result 
Re(IFT[Gf,sA(x)]) from using the spatial average, as shown in Fig. 4(c). 

The imaging system with the reference arm set in the f-f setup is able to reconstruct the 
diffraction pattern of a pure phase object, 5 both with and without the spatial averaging technique 
applied (the latter was already demonstrated in Ref. [10]). The chosen object [see Fig. 5(a)] was 
a pure phase object modulated by a Gaussian as to avoid a cosmetically disturbing DC peak 
[the analytical Fourier transform is shown in Fig. 5(b)]. An imaging scheme that is unable to 
reconstruct phase information would therefore only show the Fourier transform of the Gaussian 
which is another Gaussian. Figure 5(c) shows the result of a numerical simulation in 2D where 
the spatial averaging technique was used, and it closely follows the analytical diffraction pattern 

s Ref. [18] shows an experimental implementation of imaging the diffraction pattern of aphase object in the coinci- 
dence counting regime. 
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Fig. 5. The reconstructed diffraction pattern of a phase object with the spatial averaging 
technique, (a) shows Re(r o bj) consisting of 4 x 4 square holes (with the value T^j = —1 
inside and T \,j = +1 outside the holes) modulated by a Gaussian, (b) analytical |r o bj| 2 . (c) 
Gf.sA a f ter 2000 shots. 



of the object in Fig. 5(b) confirming that phase information is preserved. 

The image of the object can be reconstructed by merely changing the setup in the reference 
arm to the telescope setup [see Fig. 1(b)]. Figure 6 compares using (a) a pixel-like detector and 
(b) a bucket detector in the test arm for this case. In both cases the double slit is reproduced, 
albeit with the characteristic blurring of the sharp contours of the double slit apertures because 
of the finite resolution. In fact, the aperture size of the double slit (14 jj,m) is on the limit of the 
system's resolution x co h ~ 16 jim. The convergence rates as observed in the animations seem 
identical, which is confirmed by the error plot on the right [showing the image version of (12), 

i.e.e T («) = (E x |G^(x)-G T (x)| 2 ) 1/2 ]. 

From the semi-analytical calculations an interesting observation about the resolution abilities 
of the imaging scheme emerged when varying the interference filter after the PDC crystal, 
which filters a frequency range [— 5£2,5£2]. In the case considered in the analytical section 
5£2 = 0, which in practice corresponds to a < 1 nm interference filter. The numerics with 
time used <5£2 = 40£2o, i.e. a 22 nm filter. Figure 7(a) shows the case where D\ is a pixel- 
like detector; the result remains more or less unchanged as the interference filter changes and 
no noticeable difference is observed between the 22 nm and 100 nm filters. In contrast, when 
D\ is a bucket detector the result is much more sensitive; as shown in Fig. 7(b) the narrow 
interference filter gives the best resolution, while taking a broader filter gives a deteriorating 
resolution. The explanations for these different behaviours are found from considering the two 
explicit expressions for the correlations integrated over time. From Eq. (5) we namely find 
instead of Eqs. (9) and (11) 



Gt(xi,x 2 )« / d£2 



dx|r(x| -x 2 ,£2)r obj (x[)e 



-ix[-xik/f 



G T (x 2 )c< /dii /dx 1 |r(xi-x 2 ,i2)| z |r obj (x; 



coherent (pixel-like D\) (14) 



incoherent (bucket D\). (15) 



Thus, in the coherent case T(xJ — x 2 ,£2) acts at each £2 as an imaging kernel, which is shown 
in Fig. 7(d) for different values of £2. As £2 increases T(x[ — x 2 ,£2) becomes broader and 
the sidebands become pronounced. This can be understood by noting that Eq. (4) expresses 
r(x[ — x 2 ,£2) as the inverse Fourier transform of y(q,£2), which is shown in Fig. 7(c): as £2 
increases y(q,£2) becomes double-peaked, giving the sideband oscillations. All these contribu- 
tions should then be added up as the integration in £2 is carried out. It is therefore important to 
notice that r(x( — x 2 ,£2) becomes weaker as £2 is increased. Therefore, despite the fact that the 
broader sidebands in T(x[ — x 2 ,£2) imply a deteriorating resolution, the reduced peak value of 
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Fig. 6. The reconstructed image using the telescope setup with (a) D\ a pixel-like detector 
and (b) D\ a bucket detector. The black curves are numerics, the dashed red curves are 
r o bj, while the blue curves are the analytically calculated correlations in the SPWPA in 
arbitrary units [Eqs. (14) and (15)]. (c) shows ej(n) (full: bucket detector and dashed: pixel 
detector), and the red curves are fits to Eq. (13). The same parameters and notation as in 
Fig. 3. [Movie sizes: (a) 748 KB and (b) 1.25 MB]. 



r(x[ — X2,£2) implies that the end result is affected less and less as £2 is increased. Thus, the 
coherent nature of the imaging explains why the image does not change much as the bandwidth 
is increased in Fig. 7(a). In the bucket detector case the situation is completely different. Instead 
of having many contributions from different £2, the incoherence of the imaging system implies 
that the temporal integration defines a new imaging kernel as 

f d£2 , 
r B ( Xl -x 2 ) ee J ^|r(xi -x 2 ,£2)| 2 . (16) 

so Gt(x2) « / dxiTr^xi — X2,£2)|r o bj (xi) | 2 . Thus, Tb(xi — X2) uniquely determines the res- 
olution of the incoherent imaging system, and from Fig. 7(e) we see it becomes broader as 
the frequency bandwidth is increased, implying a decrease in resolution which is exactly what 
was seen in Fig. 7(b). We should stress that this result is a particular consequence of the PDC 
phase-matching conditions and does therefore not necessarily apply to ghost imaging schemes 
using other sources for the correlated beams. 

In 2D the advantage of using the bucket detector for the image reconstruction becomes more 
obvious. As an example of a reasonably complicated object, we have in Fig. 8 reconstructed a 
mask with the letters "INFM". The difference between the pixel-detector case (a) and the bucket 
detector case (b) is again related to coherent vs. incoherent imaging. 6 The semi-analytical cal- 
culations in (d) in the SPWPA confirm this; whereas the pixel-detector result is suffering from 

6 The curvature of the reproduced images is due to the Gaussian profile of the signal field impinging on the object. 
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Fig. 7. The reconstructed image from semi-analytical calculations with (a) D\ a pixel-like 
detector and (b) D\ a bucket detector, using different interferential filters (different temporal 
bandwidths). The blue curves are semi-analytical calculations in the SPWPA in arbitrary 
units [Eqs. (14) and (15)]. (c) shows \y{q x ,Q.)\ for the chosen PDC setup, (d) shows |r(x! — 
X2,£2)| for different values of SI. (e) shows Tb(xi — X2) for different interferential filters. 



a speckle-like effect [16], this effect is absent in the bucket-detector result. Note also the dif- 
ferent results obtained from a narrow and a broad bandwidth filter, as discussed in the previous 
paragraph [the numerics in Fig. 8(a)-(b) are without time, so they should be compared to < 1 
nm filter results of Fig. 8(d)]. Finally, (c) shows e(«), and it is evident that the bucket-detector 
case (full lines) converges faster than the pixel-detector case (dashed lines); from the fits we 
obtain ^ ucket /if£ ixel ~ 20. This speedup in convergence is a consequence of the bucket detector 
D\ making a spatial average over the field at the test arm detection plane, and the magnitude 
of the speedup is therefore related to the number of independent modes recorded by D\ . In ID 
this number is much smaller than in 2D which is why no significant speedup was observed in 
the ID case (Fig. 6). 

If filters are inserted in the reference arm the reconstructed object can be manipulated. As a 
simple example of this, using the telescope setup in the reference beam we have inserted a filter 
in the focal plane of the first lens [see Fig. 9(a)]. The object in the test arm was a square pattern 
?obj( x ) = [1 + cos(3xgo/2)] [1 + cos(3y^o/2)]/4 (consisting of 4 main sidebands, 2 located at 
= ±3 e -v an< ^ ^ loca ted at q/qo = ±|e v .). Filtering in the y-direction we observe only 
rolls in the x-direction [Fig. 9(b)], while removing the filter the square pattern is seen in the 
correlations [Fig. 9(c)]. This shows that the image Fourier components can be manipulated 
non-locally. 

4. Conclusion 

The ghost imaging schemes rely on two spatially correlated beams. These are generated by a 
source with a limited spatial bandwidth in Fourier space, which determines the highest Fourier 
component in the diffraction pattern of the object that can be reproduced (i.e., the imaging 
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Fig. 8. The reconstructed near field of a mask with the letters "INFM" using the telescope 
setup. In the numerical simulations of (a) D\ was a pixel-like detector, while in (b) D\ was 
a bucket detector, (c) shows the convergence rate of the two cases (same notation as in 
Fig. 6.) 10 5 shots were used, as well as a larger pump waist than usual (corresponding to 
800 ilm). (d) shows semi-analytical calculations without time (< 1 nm filter) and with time 
(22 nm filter). 




(b) 



* 




Fig. 9. Non-local filtering using the telescope setup with a bucket detector in the test arm. 
The setup (a) was as Fig. 1(b) but with a spatial filter inserted in the focal plane of the first 
lens in the reference arm. (b) shows the correlation using the filter, while (c) shows the 
correlation without filter. 15000 shots were used. 



bandwidth is limited). In turn, the resolution of the reconstructed image is limited by the char- 
acteristic near-field coherence length of the source. 

In this paper we have through theory and numerics analyzed in detail a technique (already 
presented by us in [12] for a different imaging scheme) that improves the imaging bandwidth 
by implementing a spatial average over the test arm detector position. This technique makes 
the imaging bandwidth of the reconstructed diffraction pattern virtually infinite (apart from 
limitations arising from lenses and apertures), as well as making the correlations converge much 
faster. We showed that the speed-up of this convergence is related to the number of statistically 
independent modes inside the source bandwidth. When reconstructing the object image no 
benefits to the image resolution could be made by doing a spatial average. However, by merely 
using a bucket detector collecting all photons in the test arm, the imaging becomes incoherent. 
Whether coherent or incoherent imaging is advantageous depends obviously on the problem at 



hand [16], but we showed examples where the speckle problem of coherent imaging could be 
avoided using a bucket detector. On the other hand, when using a bucket detector it turned out 
to be important to use a narrow -band interference filter of the down-converted beams, otherwise 
a strong degrading in resolution was observed. Finally, we also showed that by inserting a filter 
in the focal plane of a lens in the reference arm, the reconstructed image could be non-locally 
filtered; the filter never interferes with the field containing the object information but through 
the correlations the image is filtered nonetheless. 

We have in this paper chosen to focus on an imaging scheme using spatially entangled PDC 
beams. However, we should stress that practically all the results shown are general for any 
imaging system, quantum or classical, and are therefore relevant also for the important case 
when classical spatially correlated beams are created by splitting a thermal-like radiation on a 
beam splitter [10, 13]. In that case, e.g. the near-field correlation function (4) is instead gov- 
erned by the second-order correlation of the thermal field [10], but the overall results of this 
paper remains the same. Thus, the spatial averaging technique may still be used to extend the 
bandwidth and speed up the convergence. In this connection it should be stressed that with 
thermal-like beams the reconstructed object image has exactly the same properties as when us- 
ing PDC beams: it is coherent of one uses a pixel-like detector, while it becomes incoherent if 
a bucket detector is used. 
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